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The Hofstadter butterfly spectral patterns of lattice electrons in an external magnetic field yield some 
of the most beguiling images in physics [1-EJ. In this Letter we explore the magneto-electronic spectra of 
systems with moire spatial patterns |8j, concentrating on the case of twisted bilayer graphene. Because 
long-period spatial patterns are accurately formed at small twist angles, fractal butterfly spectra and 
associated magneto-transport and magneto-mechanical anomalies emerge at accessible magnetic field 
strengths. 

The fractal Hofstadter spectrum is a canonical example of electronic structure in a system with incommensurate 
length scales, and has fascinated physicists and mathematicians for over a half a century. The classic butterfly pattern 
is formed by the magnetic field dependent support of the eigenvalue spectrum of the Schrodinger equation for a near- 
neighbor hopping model on a square lattice. Similar but distinct [5] patterns describe the magneto-spectrum of any 
two-dimensional (2D) system of Bloch electrons. 

Quite generally magneto-Bloch Hamiltonians are block diagonalizable only when the magnetic flux through a 2D 
unit cell $ is a rational multiple of the magnetic flux quantum $o- For a = < J > o/ < I > = p/q the spectrum consists of 
continuous subbands each containing an areal density of B/qQ Q , q times smaller than the usual semiclassical Landau 
level density. Because the x and y components of cyclotron orbit centers are canonically conjugate in a magnetic 
field[10], smearing the periodic potential over the magnetic length scale I — ($o/27T-B) 1 ^ 2 , the fractal pattern of 
gaps within Landau levels becomes visible only when a is not too much larger than one. For atomic periodicity this 
condition is not met until the magnetic field strength exceeds laboratory scales by a factor of about one thousand. 
In moire systems, however, the pattern period is inversely proportional to the twist angle and can easily exceed 
I. Graphene moire systems realize Hofstadter physics at fields of a few Tesla, without recourse to the difficult and 
potentially damaging photolithographic patterning used previously ^] to realize Hofstadter physics in the lab. 

Because of the relatively weak forces between adjacent graphene layers, double layer graphene systems with a 
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variety of different stacking sequences occur in bulk graphite|12j. epitaxially grown multi-layer graphene|13). and in 
mechanically exfoliated multilayers [14!. Relative twists between layers can also be created by folding a single laver|15] 
116] . The stacking arrangement in a two-layer system can be characterized by the twist angle 9, and by a relative 
translation d. Different electronic structure aspects of the twisted bilayer have captured theoretical attention |17H23j . 
and have already spurred some experimental observations [141 124j . 

For 9 smaller than roughly ten degrees, the low energy spectrum is faithfully described by a continuum model 
obtained via an envelope function approximation [17] 123] . For small twist angles, this model shows that it is meaningful 
to describe the electronic structure in terms of Bloch bands for any 9 despite the fact that the atomic network is 
periodic only for a discrete set of angles. The Bloch bands in this description are intimately related to the moire 
pattern clearly observed in scanning tunneling microscopy measurements [13]. The moire period for bilayer graphene 
is D = a/[sin(0/2)] where a is graphene's lattice constant. Because a translation of one layer with respect to the other 
only shifts the moire pattern the electronic structure is virtually independent of d |23] except at large commensurate 
twist angles. In what follows, we therefore set d to zero. 

The continuum limit of a 7r-band tight-binding model for the twisted bilayer yields a transparent physical picture in 
which Dirac cones are coupled by a position and sub-lattice dependent interlayer hopping T(r) operator that captures 
the local coordination of the twisted honeycomb lattices. It is T(r), and not a periodic potential, which is responsible 
for the moire butterfly. Because the moire unit cell area f2 M oc 9~ 2 (see Supplementary Information), the flux through 
a moire unit cell increases rapidly as the twist angle is reduced. As in the periodic potential case, gaps open within 
Landau levels for a < 1. Because B[T] ~ 4(#°) 2 /a, significant splitting of the isolated layers Dirac Landau levels 
appear already at low magnetic fields for small 9. 

In the absence of inter-layer coupling, the spectrum consists of degenerate Dirac Landau-levels at energies iu} c s/n, 
where u c = y2u j I is the cyclotron energy, and v is the graphene sheet Dirac velocity. Because of the layer degeneracy, 
gaps between Landau levels produce quantum Hall effects at odd, rather than half odd 25 , integer filling factors. (Spin 
and valley degeneracy are left implicit throughout this article.) Coupling between the layers splits the Landau levels 
in both layers into q sub-bands and couples them together as illustrated in Fig[T]for the 9 — 2° case. It is clear that 
interlayer coupling at strong fields completely alters the spectrum. 

We now briefly derive the equations we use to evaluate the moire butterfly spectrum at rational values of a, present 
numerical results for a typical twist angle, and discuss the magneto-transport and magneto-mechanical anomalies that 



they imply. 

The low energy electronic structure of a twisted bilayer is well captured by the continuum model in which [23] 
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where <j> — 2n/3, qi = kg(0, — 1), q2 = fee(V3, l)/2, q3 = fcg(— v3, l)/2, fcg = 2fc D sin(#/2) ps fco^ with fc D being 
the Dirac momentum and w the hopping energy. Estimates based on tight-binding models for AB bilayer graphene 
suggest that w ~ HOmeV, however recent measurements suggest that w might be considerably smaller for some 
epitaxially grown layers [26 . 

In the presence of a magnetic field it is convenient to work in the Landau gauge A = B(—y,0) and express the 
Hamiltonian in the representation of the basis states \Lnay) where L = 1,2 labels the layer, n is the LL index, 
a = A, B stands for the sub-lattice, and y is the guiding center coordinate. The intra-layer part of the Hamiltonian 
is diagonal in y, however the Ti and T3 inter-layer hopping terms change y by ±A where A = \fikgl 2 /2. In the 
presence of a finite B the Hamiltonian therefore describes particles hopping on a set of one-dimensional chains. The 
Hamiltonian can be block-diagonalized in y by grouping guiding centers separated by integer multiples of A (see 
Supplementary Information). 

The guiding center chains become periodic when <5>o/& is rational, allowing a second wave vector to be introduced. 
The corresponding basis functions are constructed by writing the y-guiding coordinate as y = yo + (jnq + j)A and 
Fourier transforming with respect to m. The resulting magnetic Brillouin zone is: 



{(k u k 2 )\Q < k l =yvj? 1 < A/£ 2 ,0 <k 2 < 2vr/gA}. 



(4) 



The Hamiltonian matrix in this magnetic Bloch representation has dimension Aq times the number of Landau levels 
retained. 
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We numerically diagonalized the Hamiltonian for various twist angles accounting for inter-Landau level transitions 
(see Supplementary Information). Because a oc 9 2 /B the sub-band structure becomes more conspicuous as 9 is 
reduced. On the other hand, the band structure for very small twist angles does not have simple Dirac character even 
at B = 23\. In FigJI] we show the support of the spectrum for the intermediate case 9 = 2°. As the magnetic field 
is increased (i.e. as a is decreased) all the gaps widen. The terminology of Landau level splitting is useful as long as 
the single layer Landau levels do not overlap. Mini-gaps as large as lOmeV open up within the n = Landau level for 
B w 40T. When any one of the three tunneling processes Tj is present alone, the n = Landau level splits into two 
precisely degenerate components. The relatively large gap at the v = neutrality point which is present over a wide 
range of a in Fig[l] is a remnant of this behavior which often remains when all three hopping processes are restored. 

Since the pioneering work of Thouless et al. [5] it has been understood that the Hall conductivity <7 H (in units of 
e 2 /h) is a topological number that must be quantized when the chemical potential lies in an energy gap. Although 
the support of the spectrum as a function of field has a fractal structure, gaps in the spectra can exist continuously 
over finite ranges of field. The Landau level filling factors v at which gaps appear are characterized by two topological 
integers [5] which satisfy 

V = cr H + so. (5) 

Here 

_ il fdN\ _ n d 2 F 

A is the sample's area, fi is the area of the unit cell, N is the number of electrons in states below that gap, and T 
is the grand canonical potential. As a function of v and a, the Diophantine equation ^ has an infinite number of 
solutions: (s, er H ) = (so — m< Z: + m P) where (sq, oq) is some particular solution and m is any integer. While there is 
a simple rule to determine s is the classic Hofstadter problem[S] for the moire butterfly s and cr H must be determined 
numerically; for eaxample by plotting the energy gaps as a function of v and a. The linear dependence assured by 
equation ^ allows a straightforward identification of er H as the intercept of gap lines with the a = axis and of s as 
the gap-line slope. In Figj2]the energy gaps are plotted for 9 = 2°. 

Using Fig(2] we can identify the topological quantum numbers for every gap in Fig[TJ as illustrated for some of 
the larger gaps that appear near a m 0.3. The integers depicted in FigJT] specify the quantized Hall conductance 
in the large gaps which appear between v = 1 and v = — 1. As the electronic density is varied the quantized Hall 
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conductivity follows the non-monotonic variation +1,-1, +2, —2, +1,-1. 

As evident from equation in the case of bilayer graphene, the quantum number s can be associated with 
the chemical potential dependence of a rotational torque. Measurement of this electro-mechanical quantum number 
presents an interesting challenge to experiment. 

Our theory has intriguing consequences for magneto-transport in double-layers grown using chemical vapor depo- 
sition (CVD). This type of sample is polycrystalline in nature, characterized by graphene flakes of various sizes that 
are misoriented relative to one another. A double-layer CVD grown structure will therefore be characterized not by 
a single twist angle but by a set of 9's. In the presence of a magnetic field the Hall conductivity of each domain will 
depend on B and on the particular twist angle of the domain. Because different grains will in general have different 
Hall conductivities, chiral currents will flow along most grain boundaries. 

We note that the considerations presented here do not account for electron-electron interactions [55]. As in the 
Hofstadter butterfly, fractional quantum Hall states with fractional charge and statistics are possible (2^1 EO]- The large 
mini-gaps depicted in Figjljand the typical high mobilities of graphene multi-layers are favorable for the experimental 
observation of these fractional states in exfoliated double-layer graphene samples. 

We acknowledge a helpful conversation with Joao Lopes dos Santos, and thank Gene Mele for pointing out the 
difference between the moire pattern and the T(r) periodicity. This work was supported by Welch Foundation grant 
F1473 and by the NSF-NRI SWAN program. 
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FIG. 1: Spectrum support. The support of the spectrum as a function of a for 9 = 2° (w = HOmeV). The periodic inter- 
layer hopping amplitude results in a Hofstadter-like sub-band structure. The integers denote the Hall conductivity associated 
with the larger energy gaps between v = 1 and v — —1 for a ~ 0.3. 
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FIG. 2: Hall conductivity. A color plot for the energy gaps as a function of filling factor v and p/q for — 2° (w = HOmeV). 
The color scale corresponds to log(l + gap) in order to magnify the smaller gaps. The conspicuous straight lines satisfy the 
Diophantine gap equation |5| and determine both cr H and s. The Hall conductivity corresponds to the intercept of the line 
with the y-axis whereas s is given by its slope. 
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HAMILTONIAN 

To find the spectrum of the Hamiltonian given by equation (1) we represent it in terms of the basis states \Lnay) 
where L — 1, 2 labels the layer, n is the Landau level index, a = A, B stands for the sub-lattice and y is the guiding 
center coordinate. The intradayer part of the Hamiltonian is then 

h{9) = -uj c (er ie V^Tl\Ln + lAy) (LnBy\ + h.c.) , (7) 

Lny 

and the interdayer part is 

T= (T (0) \2nay){lmPy\ + T {R) \2nay + A)(lm(3y\ + T {L) \2nay- A){lm/3y\\ . (8) 



nmyaj: 



Here A = ^k g l 2 , 



T<°> = T\F n 



V2 



= T 2 F nm e i^(^-A) 

= T 3 F nm (^£)ei k ^+ A \ (9) 
(the Tj's are defined by equation (3) in the main text) and 

/ 777 I ,2 

F nm (z.) = {-z x +iz y ) n - m e-^C n m m (z 2 ) (10) 

for ii > m with C being the associated Laguarre polynomial. For n < m the function F can be found using 

The commensurability conditionp/g = &q/Q m B is equivalent to the condition kgA/2 = 2np/q for which the hopping 
amplitudes ^ repeat their values when y is shifted by qA. As explained in the main text, this last periodicity together 
with the 2/0 guiding center coordinate define a magnetic Brillouin zone. For each momentum k = (fci, k 2 = yo/i 2 ) in 
that zone 
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where j = 0, 1, . . . q — 1 (j is defined modulo q so that \j = q) = \ j = 0)), and 

7f> = TiF nm ( ^4) e - ik ^e-^ 



V2 



T^ R) = T,F„ m ( e i^A e ik e v 0e i^-{2j- h 
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T] L) = T 3 F nm ^)r ' '' - "\ (12i 



Because (kg£) 2 — 8irp/\/3q, the inter-layer Hamiltonian depends only on p/q. In the absence of Landau level mixing 
the splitting of each Landau level into sub-bands is therefore determined only by a. 

Inter-Landau level transitions can significantly alter the electronic spectrum[3Tl [35]. In the absence of a magnetic 
field two momentum states are effectively coupled if their energy difference is of order of E\ — max(vkg,w) or 
less. The same criteria holds also in the presence of a magnetic field. The n ~ Landau level therefore couples 
most strongly to the no « (E\/ui c ) 2 Landau level. In our calculations we retain 2no Landau levels in order to obtain 
spectra that are accurate near the Dirac point. Comparing Fig.l with results obtained neglecting Landau level mixing 
(not shown) we find that, as expected, inter-Landau level hopping is increasingly important as the magnetic field is 
reduced and as the energy is increased. 

A naive approach to implement an energy cutoff would be to keep all states whose Landau level index is less than 
2uq. The problem with this approach is that it introduces a fake zero energy state irrespective of the size of no. We 
therefore make the energy cutoff such that only one sub-lattice of the no Landau level is retained. This approach 
shifts the numerical error to high energies, of the order of ^/riouj c . 



MOIRE UNIT CELL 

The inter-layer hopping in twisted double layer graphene is akin in some respects to a periodic potential. In a real 
space representation it is described by the 4x4 matrix (layer x sub-lattice) 



' T(r) ^ 



(13) 



yTt(r) j 

where T(r) is given by equation (2) in the main text. The AA entry of the T-matrix is depicted in Figj3] (similar 
figures are obtained for the other entries of T). The direction and size of an arrow at point r correspond, respectively 
to the phase of T AA (r) and to its magnitude. The red dots mark the lattice associated with the moire pattern whereas 
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the green circles mark the lattice induced by T(r). Interestingly, the moire unit cell area 



ft M = I6n 2 /V3k 



(14) 



is six times larger than the area of the moire pattern unit cell[ 
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FIG. 3: Moire unit cell. The AA entry of the hopping amplitude is plotted as a function of position (in units of fc^" 1 ). The 
space dependent eigenvalues of H T lead to the moire pattern period marked by the red dots, however the spatial variation of 
the phase of T results in a larger period denoted by the green circles. It is this latter periodicity that determines the moire 
unit cell area Q M . 



